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Collective modes of a spin-orbit-coupled superfluid Fermi gas in a two-dimensional 
optical lattice: a comparison between the Gaussian approximation and the 

Bethe-Salpeter approach 

Zlatko Koinov, Rafael Mendoza 
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A functional integral technique and a Legendre transform are used to give a systematic deriva¬ 
tion of the Schwinger-Dyson equations for the generalized single-particle Green’s function and the 
Bethe-Salpeter equation for the two-particle Green’s function and the associated collective modes of 
a population-imbalanced spin-orbit-coupled atomic Fermi gas loaded in a two-dimensional optical 
lattice at zero temperature. The collective-mode excitation energy is calculated within the Gaussian 
approximation, and from the Bethe-Salpeter equation in the generalized random phase approxima¬ 
tion assuming the existence of a Sarma superfluid state. It is found that the Gaussian approximation 
overestimates the speed of sound of the Goldstone mode. More interestingly, the Gaussian approx¬ 
imation fails to reproduced the rotonlike structure of the collective-mode dispersion which appears 
after the linear part of the dispersion in the Bethe-Salpeter approach. 

We use the Gaussian approximation and the Bethe-Salpeter approach to investigate the speed 
of sound of a balanced spin-orbit-coupled atomic Fermi gas near the boundary of the topological 
phase transition driven by an out-of-plane Zeeman field. It is shown that within both approaches, 
the minimum of the speed of sound is located at the topological phase transition boundary, and this 
fact can be used to confirm the existence of a topological phase transition. 

PACS numbers: 03.75.Kk, 03.75.Ss, 67.85.Lm 


I. INTRODUCTION 

Topological superfluidity is an interesting state of mat¬ 
ter, partly because it is associated with quasiparticle exci¬ 
tations which are Majorana fermions X The basic physics 
behind the emergence of the Majorana fermion excita¬ 
tions is the existence of s-wave superfluidity, nonvanish¬ 
ing spin-orbit coupling (SOC), and Zeeman splitting. In 
this context, calculating the dispersion of the collective 
modes when the pseudospin of atoms can couple with 
not only the effective Zeeman field, but also with the or¬ 
bital degrees of freedom of atoms is an important and 
interesting problem by itself. 

The recent experimental breakthroughs in realization 
of spin-orbit coupling (SOC)2r— in ultracold atomic 
gases have opened up possibilities for investigating the 
phase diagram, the single-particle and the collective- 
mode excitations of superfluid Fermi gases in optical 
lattices in the presence of SOC and Zeeman fields. It 
is widely accepted among the optical-lattice community 
that the attractive Fermi-Hubbard model captures the 
s-wave superfluidity of cold fermion atoms in optical lat¬ 
tices. According to this model, two fermion atoms of 
opposite pseudospins on the same site have an attractive 
interaction energy U, while the probability to tunnel to 
a neighboring site is given by the hopping parameter J: 

Hjj — ^ ^ ^ ^ ^ /GTRer- 

<i,j>,cr i i,a 

Here, J a is the tunneling strength of the atoms between 
nearest-neighbor sites, and fij jCr = ip\ is the density 
operator on site i. The Fermi operator ip\ a (^>j jCr ) creates 
(destroys) a fermion on the lattice site i with pseudospin 


projection a =t,4-- The symbol '}2 <ij> means sum over 
nearest-neighbor sites of the two-dimensional (2D) lat¬ 
tice. 

In this paper we assume the existence of nonvanishing 
Rashba SOC in the xy plane and a Zeeman field along 
the z direction, so the Hamiltonian of the system is H = 
Hu + Hsoc + Hz- In the case of a 2D optical lattice the 
SOC part of the Hamiltonian is given b y 13 i 14 


Hsoc = -tX X ^ x d o) z ( 

<i,j> ' 


tpiA 


where A is the Rashba SO coupling coefficient, = 
(a x , & y , Cz), a x,y,z are the Pauli matrices, and d,;j is a 
unit vector along the line that connects site j to i. 

The out-of-plane Zeeman field is described by the term 
H z : 


H z 



The simplest way to study the collective excitations of 
the above-mentioned model is to apply functional integral 
technique which requires the representation of the Hub¬ 
bard interaction in terms of squares of one-body charge 
and spin operators. It is possible to resolve the Hub¬ 
bard interaction into quadratic forms of spin and electron 
number operators in an infinite number of ways by ap¬ 
plying the Hubbard-Stratonovich transformation. If no 
approximations were made in evaluating the functional 
integrals, it would no matter which of the ways is cho¬ 
sen. When approximations are taken, the final result 
depends on a particular form chosen. 
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One of the most common ways to apply the Hubbard- 
Stratonovich transformation is to introduce the energy 
gap as an order parameter field, which allows us to inte¬ 
grate out the fermion fields and to arrive at an effective 
action. The next steps are to consider the state which 
corresponds to the saddle point of the effective action, 
and to write the effective action as a series in powers of 
the fluctuations and their derivatives. The exact result 
can be obtained by explicitly calculating the terms up 
to second order in the fluctuations and their derivatives. 
This approximation is called the Gaussian approxima¬ 
tion. 

In the case of vanishing SOC and h~ = 0, the Gaus¬ 
sian approximation provides the following equation for 
the collective mode dispersion w(Q) of homogeneous su¬ 
perfluid Fermi gases at a zero temperature^ 

0 = 1 + 1/ (Xy i7 + A,z) + U 2 — Jyj) > (1) 

where / 7i7 = / 7i7 (w, Q), A,; = /;,( ui,Q) and J 7i j = 
J 7i ;(w, Q) are defined in the Appendix. 

Instead of introducing the energy gap as an order pa¬ 
rameter field, one can transform the quartic terms to a 
quadratic forms by introducing a four-component boson 
field which mediates the interaction of fermions. This 
approach is similar to the situation in quantum electro¬ 
dynamics, where the photons mediate the interaction of 
electric charges. This similarity allows us to apply the 
powerful method of Legendre transforms, to derive the 
Schwinger-Dyson±&iI (SD) equation G" 1 = G^" 1 - E, 
and the Bethe-Salpeter— (BS) equation [AT®” 1 — /]<P = 
0 for the poles of the single-particle Green’s function, G, 
and the poles of the two-particle Green’s function, AT, 
respectively. Here, G+) is the free single-particle propa¬ 
gator, E is the fermion self-energy, I is the BS kernel, and 
the two-particle free propagator A++ = GG is a product 
of two fully dressed single-particle Green’s functions. The 
kernel of the BS equation is defined as a sum of the direct 
interaction, Az = ST, F /SG, and the exchange interaction 
I exc = ST, H /SG, where E^ and T, H are the Fock and 
the Hartree parts of the fermion self-energy E. Since the 
fermion self-energy depends on the two-particle Green’s 
function, the positions of both poles must to be obtained 
by solving the SD and BS equations self-consistently. 

It is widely accepted that the generalized random 
phase approximation (GRPA) is a good approximation 
for the collective excitations in a weak-coupling regime, 
and therefore, it can be used to separate the solutions of 
the SD and the BS equations. In this approximation, the 
single-particle excitations are obtained in the mean field 
approximation (or by solving the Bogoliubov-de Gennes 
equations in a self-consistent fashion); while the collec¬ 
tive modes are obtained by solving the BS equation in 
which the single-particle Green’s functions are calculated 
in Hartree-Fock approximation, and the BS kernel is ob¬ 
tained by summing ladder and bubble diagrams. The 
BS equation in the GRPA has been used to obtain the 
collective-mode spectrum of an imbalanced Fermi gas in 
a deep optical lattice *2r— 


In the GRPA, the equation (JT]) for the collective mode 
dispersion w(Q) of homogeneous superfluid Fermi gases 
at a zero temperature is replaced by the following one >2^ 


0 — 1 + (A,z + / 7 , 7 + An,m) U + (A,zA 7>7 — i 
A, m + A,/Ara,m '/y, rn + An,mAy )7 )t/ + ( 

+ 2A, m A 7j /A 7 , m A,iJ 7 m / 7l7 A m T Il,lIm,mI'y,~t)U . 

( 2 ) 


In Fig. Q] and Fig. [2J we present the collective¬ 
mode dispersions ui(Q x ) in a one-dimensional (ID) op¬ 
tical lattice, and oj(Q x ,Q y ) in a 2D lattice along the 
Qx = Qy = Q direction, both obtained by numerically 
solving Eqs. m and ©• The speed of sound is defined by 
the slope of the Goldstone-mode dispersion in the limit 
Q —> 0. When the Gaussian approach is used, the speed 
of sound in the ID case is u = 1.30 Ja/h 1 while the GRPA 
provides u = 0.96 Ja/h. In the 2D case, the correspond¬ 
ing speeds are u = 2.00 Ja/h, and u = 1.30 Ja/h, respec¬ 
tively. In both cases, the speed of sound is overestimated 
by the Gaussian approximation. In a diagrammatic lan¬ 
guage, Eq. © can be derived by summation of the infi¬ 
nite sequences of graphs in the ladder approximation! 24 
Thus, one can say that the speed of sound decreases 
when the exchange interaction, represented by bubble 
diagrams, is taken into account. Another interesting fact 
is that in ID case both approximations provide the ex¬ 
istence of a roton minimum, while the inset (b) in the 
Fig. [2] shows no roton minimum within the Gaussian 
approximation. 

In what follows, we study the collective-mode disper¬ 
sion of species of Fermi atoms with equal, or imbalance, 
population in two pseudospin states loaded in a 2D op¬ 
tical lattice in the presence of both Zeeman field and 
nonvanishing Rashba type of spin-orbit coupling. To the 
best of our knowledge, the Gaussian approximation is 
the only approximation that has been used to obtain the 
speed of sound in the presence of the SOC and the Zee- 
man field effects^ - — In a view of the facts that: (i) in a 
lattice system the bubble diagrams induce the instability 
due to the charge-density-wave (CDW) fluctuations, (ii) 
because of the strong CDW fluctuations the collective¬ 
mode spectrum has a characteristic roton-like structure 
which lies below the particle-hole continuum (for a de¬ 
tailed discussion of the effects of the CDW fluctuations 
on the stability of the collective modes, see Ref. [33l]+ and 
(iii) the speed of sound is overestimated in the Gaussian 
approximation, one may well ask what would be the dif¬ 
ference between the collective-mode dispersion and the 
corresponding speed of sound in the presence of the SOC 
and Zeeman fields when calculated within the Gaussian 
approximation, and from the BS formalism by summing 
infinite sequences of ladder and bubble diagrams. 

The main goal of the present study is to obtain the 
BS equation for the collective modes in the GRPA which 
takes into account the SOC and the Zeeman effects. In 
Sec. II, we derive the BS equation for the poles of the 
two-particle Green’s function, which allows us to ob- 
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FIG. 1: Collective mode dispersions of the Goldstone mode 
in a one-dimensional optical lattice in the weak coupling limit 
obtained by applying the Gaussian approximation (triangles), 
and by the GRPA (circles). The lattice constant is a. We set 
the filling factors to be ,ft = ft = 0.25, and the strength of 
the interaction is U = 2 J. The mean field superfluid gap and 
the chemical potential are A = 0.409J and n = 0.624J. 



Q (units of n/a) 


FIG. 2: Collective mode dispersions lo(Q x ,Q v ) of the Gold- 
stone mode in a two-dimensional optical lattice along the 
Q x = Q y = Q direction in the weak coupling limit obtained by 
applying the Gaussian approximation (triangles), and by the 
GRPA (circles). The rotonlike structure is seen in the GRPA 
spectrum (inset (a)), while the roton minimum does not ap¬ 
pear in the Gaussian approximation (inset (b)). The filling 
factors, and the strength of the interaction are ft = fl = 0.25 
and U = 4.5 J. The mean field superfluid gap and the chemi¬ 
cal potential are A = 1.334J and /x = 2.234J. 


tain numerically the dispersion of the collective modes 
in the presence of both the Zeeman field and the Rashba 
types of SOC. In Sec. Ill, the roton-like structure of 
the collective-mode spectrum of an imbalanced 2D Fermi 
gas in the presence of SOC is obtained by solving the 
BS equation in the GRPA. It turns out that the Gaus¬ 
sian approximation and the BS approach, both provide 
a very similar slope of the linear part of the Goldstone 
dispersion, but there is no roton minimum within the 
Gaussian approximation. In the same Section, we study 
the speed of sound as a function of the strength of the 
Zeeman field near the topological quantum phase tran¬ 
sition boundary from a nontopological superfluid state 
with a fully gapped fermionic spectrum to a topological 
superfluid state. 


II. THE BETHE-SALPETER EQUATION 

A. The functional-integral formulation of the 
Hubbard model 

We consider an imbalanced mixture of an atomic Fermi 
gas of two hyperfine states (described by pseudospins 
a =t, |) with contact interaction loaded into a 2D square 
optical lattice. In the imbalanced case, there are differ¬ 
ent amounts of Mf and Mj. atoms in each state achieved 
by considering different chemical potentials ^t an d Ml- 
The total number of atoms M = + Mj. is dis¬ 

tributed along N sites, and the corresponding filling fac¬ 
tors = Mt,i/N are smaller than unity. 

From a theoretical point of view, the corresponding 
expressions for the Green’s functions cannot be evalu¬ 
ated exactly because the interaction part of the Hubbard 
Hamiltonian is quartic in the fermion fields. The simplest 
way to solve this problem is to apply the so-called mean- 
field decoupling of the quartic interaction ; 13 i 14 To go be¬ 
yond the mean-field approximation, we apply the idea 
that one can transform the quartic term into quadratic 
form by making the Hubbard-Stratonovich transforma¬ 
tion for the fermion operators. In contrast to the previ¬ 
ous approaches, wherein after performing the Hubbard- 
Stratonovich transformation the fermion degrees of free¬ 
dom are integrated out; we decouple the quartic problem 
by introducing a model system which consists of a multi- 
component boson field A a interacting with fermion fields 
tfj' and if. 

The functional-integral formulation of the Hubbard 
model requires the representation of the Hubbard inter¬ 
action Hu in terms of squares of one-body charge and 
spin operators. It is known that when approximations 
are made, the final result depends on the particular form 
chosen. Thus, one should check that the results obtained 
with a certain Hubbard-Stratonovich transformation are 
consistent with the results obtained with canonical mean- 
field approximation. It can be seen that our approach 
to the Hubbard-Stratonovich transformation provides re¬ 
sults consistent with the results obtained with the mean- 
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field approximation, i.e. one can derive the mean-field 
gap equation using the collective-mode dispersion uj(Q) 
in the limit Q —> 0 and ui —> 0. 

The Green’s functions in the functional-integral ap¬ 
proach are defined by means of the so-called generat¬ 
ing functional with sources for the boson and fermion 
fields, but the corresponding functional integrals can¬ 
not be evaluated exactly because the interaction part of 
the Hubbard Hamiltonian is quartic in the Grassmann 
fermion fields. However, we can transform the quartic 
terms to a quadratic form by introducing a model sys¬ 
tem which consists of a four-component boson field A a (z) 
(a = 1,2,3,4, z = (!■;,«), 0 < v < (3 = (fc B T) _1 , 

ks is the Boltzmann constant) interacting with fermion 

fields = 4 A(y)/\/2 and t/>(x) = ^>(x)/y/2, where 


'h(x) 


/ 4>t( x ) \ 
i>d x ) 

ipUx) ’ 

V / 


&(y) = (V'](2/)V , |(y)V’t(y)V'4(2/)), 


( 3 ) 

The field operators ([3]) allow us to define the gener¬ 
alized single-particle Green’s function by using a ten¬ 
sor product of these two matrices. The corresponding 
Green’s function, represented by a 4 x 4 matrix, includes 
all possible thermodynamic averages: 


G(x 1 ;y 2 ) = - < T u ($(xi) <g> ^( 2 / 2 )) > • (4) 

I 


The action of the above-mentioned model system is as- 

( F''\ ( R'l 

sumed to be of the following form S = Sq + Sq + 
S^ F ~ B \ where Sq F) = i/j(y)G ( ' 0 ' > ~ 1 (y : x)tp{x), Sq B ^ = 
\A a {z)D { d ) ~ 1 (z,z')A p {z'), S= ^{y)f ( a\y,x \ 

z)i/j(x)A a (z). The action S^" 1 describes the fermion part 
of the system. The generalized inverse Green’s function 
of free fermions G*- 0 ' 1 (y; x) is given by the following 4x4 
matrix: 


G^~ 1 (y,x) = 

^2 exp [*k.(ri - iv) - ui m (u - it')] G^}~2 k, uv m ), 

k,o; m 


where G^ 1 (k,uo m ) = — G^ 1 (—k and 

-G^ )_ 1 (-k, —iu m ) = G{ 0 ^ 1 {k 7 iui rn ). The symbol 

is used t° denote j3~ l Ym (f° r f erm i° n fields 
uj m = (27r//3)(m + 1/2); m = 0, ±1, ±2,...). In the case 
of the population-imbalanced Fermi gas with a Rashba 
SO coupling and an out-of-plane Zeeman field, the non¬ 
interacting Green’s function is<2£ 


G^-\k,zoj m ) 


( iuj m — £f(k) — h z 2 A (sin k x + 1 sin k y ) 0 

2 A (sin k x — 1 sin k y ) iuj m — £j.(k) + h z 0 

0 0 iu) m + £ t (k) + h z 

\ 0 0 2 X(smk x + ismky) 


. » . ) 

— 2 A (sin k x — 1 sin k y ) 1 

lUm +Ct(k) - h z ) 


( 5 ) 


Here, Ct,l(k) = 2 4,(1 — cos k x ) + 2 J-|-,l(l — cos k y ) — P77 is the tight binding form of the electron energy (the lattice 

constant a = 1 ). 

/ D\ 

The action Sq ; describes the boson field which mediates the fermion-fermion on-site interaction in the Hubbard 

/ D\ 

Hamiltonian. The bare boson propagator in Sq ; is defined as: 


D^°\z, z') = 6(v 


v')USj,j' 


/ 0 1 0 0 \ 
1 0 0 0 1 
0 0 0 0 
V 0 0 0 0 / 


The Fourier transform of this boson propagator is given by 


D (0 \z,z') = l^^e{ , [ k -( r ^- r i'H' 1 ( W )]}S ( 0 ) (k),S ( 0 ) (k) = 


k ojp 


/ 0 

u 

0 

0 

u 

0 

0 

0 

0 

0 

0 

0 

V 0 

0 

0 

0 


( 6 ) 


The interaction between the fermion and the boson fields is described by the action S^ F B \ The bare vertex 
Va\yi’,X 2 I z) = T ( 'a\i\,ui\i 2 1 U 2 \ j,v) = d(ui - u 2 )5(ui - tO^ui^iuF 1 ^(a) is a 4 x 4 matrix, where 

f ( 0 ) (a) = -(70 + a z )S a 1 + ^(70 - a z )S a2 + ^(a x +ia y )6 a 3 + ^(ce x - ia y )5 a 4- 


( 7 ) 
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The Dirac matrix 70 and the matrices on are defined as (when a four-dimensional space is used, the electron spin 
operators cq has to be replaced by ): 



/I 

0 

0 

0 N 





0 

1 

0 

0 


O 

b 


7o = 

0 

0 

-1 

0 

, OLi = 

^ 0 CTyCTiCTy J 

,i = x,y,z 


\0 

0 

0 

- 1 ) 




The relation between the Hubbard model and our model system can be demonstrated by applying the Hubbard- 
Stratonovich transformation for the fermion operators: 


Dy[A\ exp \ip(y)f^\y,x\z)ip(x)A a (z) 


= exp 




The functional measure Dy[A\ is chosen to be: 

Dfi[A] = DAe-^ A ^ z)D< '°y^’ z ' )A ^ z '\ j [m[A] = 1 . 

According to the field-theoretical approach, the expectation value of a general operator 0(u ) can be expressed as 
a functional integral over the boson field A and the Grassmann fermion fields ip and ip: 

1 


< T u (0(u)) >= 


Dy[ip,ip,A]0(u)ex p J a (z)A a (z) - M(ip,ip) 


|J=M=(b 


Z[J,M] 

where the symbol < ... > means that the thermodynamic average is made. The functional Z[J, M] is defined by 


( 8 ) 


Z[J,M] = / Dy[ip,ip,A]ex p J a (z)A a (z) - M(ip,ip) 


(9) 


where the functional measure Dfj,[ip, ip, A] = DADipDip exp ( S ) satisfies the condition j D/j,[ip, ip, A] = 1. The quantity 

J a {z) is the source of the boson field, while the sources Mij(y;x) of the fermion fields are included in the M{ip,ip) 
term : 

M{ip,ip) =ip\(y)M n (y,x)ip t {x)+ipl(y)M 2 i{y;x)ip t (x)+ip\(y)M 12 (y;x)ipi{x) 

+ipl(y)M 2 2 {y; x)ipi(x) -t -ip^(y)M 31 (y-x)ip t (x) + ip i {y)M 41 (y, x)ip t {x) 

-M Pt(y)M 32 (y; x)ipi{x) + ipi(y)M 42 (y-x)ip±{x) 

+ip\(y)M 13 (y; x)x)ip\(x) + ip\{y)M 23 (y, x)ip\(x) + ip\(y)M u (y, x)ip\{x) 

+ipl(y)M 2i (y; x)ip\{x) + ip^(y)M 33 (y; x)ip\(x) + ipi{y)M 43 (y; x)ip\(x) 

+ip t (y)M 34 (y-x)ip\{x) + ipi(y)M 44 (y; x)ip\{x). (10) 

We shall now use a functional derivative 8/8M{ 2; 1) = 8/8M n2 , ni (y 2 ;x 1 ), where 1 = {m,x 4 } and 2 = {n 2 ,y 2 } 
are complex indices; depending on the spin degrees of freedom, there are sixteen possible derivatives. By means of 
the definition ©, all Green’s functions can be expressed in terms of the functional derivatives with respect to the 
corresponding sources of the generating functional of the connected Green’s functions W[J,M] = In Z[J,M], Thus, 
we define the following Green’s and vertex functions which will be used to analyze the collective modes of our model: 
The Boson Green’s function is D a p(z, z') is a 4 x 4 matrix defined as 

S 2 W 

D a p(z, z ) = ~ j , , ,■ 

8J a (z)8Jp(z') 

The generalized single-fermion Green’s function G nin2 (x 1 ; y 2 ) is the matrix defined by Eq. (U) whose elements are 

G nin2 (xi, t/ 2 ) — 8W/8M n2ni (y 2 , x 4 ). 

Depending on the two spin degrees of freedom, f and 4-, there exist eight ’’normal” Green’s functions and eight 
’’anomalous” Green’s functions. We introduce Fourier transforms of the ’’normal” G aii(T2 (k, u 4 — u 2 ) = < 
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Tu(t/Vi, k(ui)Vd 2 ,k( u 2 )) >, and ’’anomalous” F ait(T 2 (k, u\-u 2 ) = - < T„(Vv llk (Mi)Vv 2 ,-k(« 2 )) > one-particle Green’s 
functions, where { 04 , < 72 } =t,4-- Here and ip^ k {u), V4,k( w ) are the creation-annihilation Heisenberg 

operators. The Fourier transform of the generalized single-particle Green’s function is given by 


G(l;2) 


Tr Ek E^ m ex P(* [ k - ( r u - D 2 ) - WmOi - « 2 )]} 


G(k, UV m ) F{k,lOJm ) \ 
F^(k,iu>m) -G(~k, -uo m ) ) 


( 11 ) 


Here, G and F are 2x2 matrices whose elements are G ai ^ 2 and F ait „ 2 , respectively. 


The two-particle Green’s function K 


n i,xi n 3 ,y 3 

n 2 ,j / 2 ri 4 ,x 4 


K 


ni.xi n 3 ,y 3 
n 2 l y 2 ri 4 ,X 4 


= K 


1 3 

2 4 


is defined as 
S 2 W 


_ 8G nin2 (zi; 1 / 2 ) 

SM n 2 ni (y 2 ', Xi)SM n3ni (y 3 ; x 4 ) 8 M n 3 Tli (y 3 \ x 4 )' 

The vertex function T a (2; 1 | z) for a given a is a 4 x 4 matrix whose elements are: 

8G-} ni {i 2 ,u 2 ]ix,U4) 


r„(*2,u 2 ;*i,ui | v , j) 

r 


8Mz') 


D n a ( z 'i z )- 


Next, we shall obtain the corresponding equations of the boson and fermion Green’s functions. The poles of these 
Green’s functions provide the single-particle and the two-particle excitations. 

It is known that the fermion self-energy (fermion mass operator) E(l; 2) can be defined by means of the SD equations. 
They can be derived using the fact that the measure Dij[ijj, %/j, A] is invariant under the translations ip —>■ ip + Si/j and 
A ->■ A + 5A: 


D^-\z,z')Rp{z') + \jr (G(l;2)f(°>(2; 1 | *)) + J a (z) = 0, 


( 12 ) 


G _1 (l; 2) - G ( °J _1 (1; 2) + E(l; 2) + M( 1; 2) = 0, 


(13) 


where R a (z ) = SW/8J a (z) is the average boson field. The fermion self-energy E, is a 4 x 4 matrix which can be 
written as a sum of Hartree T, H and Fock H F parts. The Hartree part is a diagonal matrix whose elements are: 

E ff (ii,ui;* 2 ,-u 2 )„ 1 n 2 = |f (ii, ui; * 2 , u 2 \j, v) ni „ 2 (j,v, j’, v') 

r^ 0) (*3,M3;*4,M4|/ ) V ) 713714 Gn 4 7i 3 (u,U4;i3,U3). (14) 


The Fock part of the fermion self-energy is given by: 


S F (ii,u 1 ;i 2 ,u 2 ) niri2 = —ri 0) (ii,ui;* 6 , u 6 jj,v) 


5 D a0(j: v 'if> v ') 


p( 0 ) 

1 0 


( 24 , U 4 ', 25 , U 3 \j ,27 ) n4U5 K 


n 5 ,i 5 ,u 5 n 3 ,i 3 ,u 3 


n-i 

^n 3 n 2 


(* 3 > u 3 ; i 2 , u 2 ). 


(15) 


The Fock part of the fermion self-energy depends on the two-particle Green’s function K ; therefore the SD equations 
and the BS equation for K have to be solved self-consistently. 

Our approach to the Hubbard model allows us to obtain exact equations of the Green’s functions by using the field- 
theoretical technique. We now wish to return to our statement that the Green’s functions are the thermodynamic 
average of the T u -ordered products of field operators. The standard procedure for calculating the Green’s functions, is 
to apply the Wick’s theorem. This enables us to evaluate the T u -ordered products of field operators as a perturbation 
expansion involving only wholly contracted field operators. These expansions can be summed formally to yield different 
equations of Green’s functions. The main disadvantage of this procedure is that the validity of the equations must 
be verified diagram by diagram. For this reason we shall use the method of Legendre transforms of the generating 
functional for connected Green’s functions. By applying the same steps as in Ref. [37}, we obtain the BS equation of 
the two-particle Green’s function, the Dyson equation of the boson Green’s function, and the vertex equation: 


K~i ( n 2 , i 2 , U 2 713,13,163 

l 77 - 1 ,21,221 224,24,224 


^(0) —1 ( n 2,^2, U 2 n 3 ,i 3 ,U 3 
\ ni, 21 , 2/1 224 , 24,224 


J { 742, *2, U 2 n 3 , i 3 , U 3 

yni : iii u i 714,14,164 


(16) 
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D n 


p{z, z') = 0, z') + (z , z")ILy S (z'‘ 


')Dfi 


(z,z% 


(17) 


p /• | \ _p(0)/- • | \ , T n 2,l2,V>2 71 3 ,! 3 ,U 3 \ 

T a (t 2 ,U 2 ,H,U 1 \z) n2 n 1 -T a (t 2 ,«2,*l,«l I *)n a m + / ( ru,n,U A ) 


K (°) 


n3) l 3, ^3 TIq,Ig,Uq 
^■4? ^4? ^4 ^5} ^ 5 , ^5 


r ck(^ 6? ^6? ^5? ^5 | '2’)r 


Here, 


ir< 0 > 


712, *2, M2 M 3 ,l 3 ,M 3 
711,21,1^1 774 , 24,724 


= G n2rl3 (?2, M3! *2, M2)G ri4ni (*4, M4; *1, Ml) 


is the two-particle free propagator constructed from a pair of fully dressed generalized single-particle Green’s functions. 
The kernel I = 6T,/6G of the BS equation can be expressed as a functional derivative of the fermion self-energy £. 
Since £ = T, H + £ F , the BS kernel I = I exc + Id is a sum of functional derivatives of the Hartree Y> H and Fock S F 
contributions to the self-energy: 


Ie 


7l2,«2,M2 713,i 3 , M3 \ <5XT(l2, 7X25*1, Ml) 


77-1,21,721 774 , 24,724 


SG n3 

n 4 (*3, M3; *4, M4) 


Id. 


712, *2, M2 713,13, M3 
771, Zl, 721 774, 24, 724 


(5E f (i 2 ,m 2 ;ii,mi) 

^ 2 ni 

^G n3 n A (23, 723, 24, 724) 


The general response function n in the Dyson equation m is defined as 


,K 


II a/ 3 (z,z') = f^ 0) (li,Mi;i 2 ,M 2 | z) 

The functions D, K and T are related by the identity: 


7l2,l2,M 2 ti 3 , i 3 , M 3 (p(o)^- 
771, Z1,721 774, 24, 724 


T(^3 5^35^45^4 | % )n,3n4' 


(18) 


(19) 


7 ^( 0 ) ( ^ 25 ^ 25^2 ^ 35 ^ 35^3 \ p / • * I /\ 7 ~) / / \ 

^ ( m,*i,Mi 714,1 4 , M 4 ) r /3(74,M4,l 3 ,M 3 \z )n 4 n 3 Dp a (z , z) 


= K 


712,12, M2 71 3 ,13, M 3 
7715 Z15 721 774, 24, 724 

3-1/ 


r^ 0 ) (i 4 ,M 4 ;i 3 ,M 3 | 


( 0 ), / 


( 20 ) 


By introducing the boson proper self-energy P a J (z, z') = 11^(2, z')+D ( ^ (z. z'). one can rewrite the Dyson equation 
m for the boson Green’s function as: 


D~^z, z') = D^~\z, z 1 ) - P a p{z, z 1 ). 

The proper self-energy and the vertex function T are related by the following equation: 

P a p(z,z') = ^Tr (y 1 ,x 2 \z)G(x 2 ,y3)I'p{y3,X4\z')G(x4,yi) 

= hi°)(n , Mi,l2,M2 | z)n in2 G n2 n 3 (l2, M2, I3, M3)T ( g(l 3 , M 3 , I4, M4 | Z )n 3 ri 4 G n ^ ni (14, M4, i \, Ml). 

It is also possible to express the proper self-energy in terms of the two-particle Green’s function K which satisfies 
the BS equation K~ l = — Id , but its kernel Id = <5E F / 5G includes only diagrams that represent the direct 

interactions: 


( 21 ) 


P a p(z,z') = r®(li,Mi;i 2 ,M 2 | Z) ni n 2 l 

= f( 0 ) ( n )K ( n2,r P V n Zi r 3'’ V ' \ f( 0 ) to\ 

rai712 ^ ' yni,rj,v 714, ry, v 1 ) n 3 n ^’' 


712,12, M2 m 3 ,13 , m 3 
71i,ii,Ml 7l4,l4,M4 


T^ 0) (i3,M 3 ;i4,M4 | Z')n 3 n 4 


One can obtain the spectrum of the collective excitations as poles of the boson Green’s function by solving the Dyson 
equation (ED. but first we have to deal with the BS equation for the function K. In other words, this method for 
obtaining the collective modes requires two steps. For this reason, it is easy to obtain the collective modes by locating 
the poles of the two-particle Green’s function K using the solutions of the corresponding BS equation. 
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As we have already mentioned, the BS equation and the SD equations have to be solved self-consistently. In what 
follows, we use an approximation which allows us to decouple the above-mentioned equations and to obtain a linearized 
integral equation for the Fock term. To apply this approximation, we first use Eq. to rewrite the Fock term as 

£ F (*i,ui;t 2 , ^ 2 ) niri2 -rL 0) (ir , Wi; * 3 ) u 3 \j, v) nin3 D a p(j, v; j\ v')Gn 3ni (i 3 , u 3 ; i 4 , U 4 )Tp(i 4 , u 4 ; i 2 , u 2 \j' 5 ^ )n4?l2 5 (22) 

and after that, we replace D and T in (l22l) by the free boson propagator D and by the bare vertex I^ 0 ', respectively. 
In this approximation the Fock term assumes the form: 


£0 (* 1 ) Ui; * 2 , U2)nin a = -ri 0 ) (*i, Ui; i 3 , u 3 \j, v) nin3 D { °}(j , v; j', u')f i 0) (z 4 , u 4 ; i 2 , u 2 |/, « , ))n 4 n 2 G ! rt 3„4(i3, u 3 ; l 4 , u 4 ) = 


- U6 ilii2 6(ui - u 2 ) 


a/3\ 

( 0 Gi 2 (1;2) 0 -G 14 (l;2) 

G 2 i(1; 2 ) 0 —G 23 (l; 2) 0 

0 —G 32 (l; 2) 0 G 34 (l; 2) 

V —G 41 (l;2) 0 G 43 (l; 2) 0 


(23) 


The total self-energy is £(zi, u\\ 22 , U 2 ) = u\\i^ U 2 ) + XT (ii, 12 , ^ 2 ), where 


T, H (ii,ui m ,i 2 ,u 2 ) — — (Jq^^iti — u 2 ) 


(G 22 - G 44 0 0 0 

0 Gn — G 33 0 0 

0 0 G 44 - G 22 0 

Vo 0 0 G 33 - G 11 


(24) 


r 


where Gq = Gq(l;2) = Gij(ii,ui;i 2 ,u 2 ). 

The contributions to £(q, u±; i 2 , u 2 ) due to the ele¬ 
ments on the major diagonal of the above matrices can 
be included into the chemical potential. To obtain an 
analytical expression for the generalized single-particle 
Green’s function, we assume two more approximations. 
First, since the experimentally relevant magnetic fields 
are not strong enough to cause spin flips, we shall assume 
Gi 2 = G 2 i = G 34 = G 43 = 0. Second, we neglect the fre¬ 
quency dependence of the Fourier transform of the Fock 
part of the fermion self-energy. Thus, the Dyson equa¬ 
tion for the generalized single-particle Green’s function 
becomes: 


G~ 



/ G^' 

-1 

r (0)-l 

(_T 12 

0 

-Aq 



‘(i;2) = 

/"i( o)- 
( “ r 21 

0 

-1 

r (o)-i 

^22 

^ii ,12 

Aq ; q 

°33 

0 

r (o)~ 

(j 31 

1 

? 


V — Aq 

5^2 

0 

r-(o)-i 

ljr 41 

sy(0)~ 

Lt 44 

1 J 


G<”>- = 

r (0)-l 

u a 






(25) 

(i; 

2). In the population-balanced 

= 

A 0 <5(r 

h 

- d 2 ), 

where 

A 0 is 

a 

con- 


stant in space. Physically, it describes a superfluid 
state of Cooper pairs with zero momentum. Superfluid 
states of Cooper pairs with nonzero momentum occur 


J 


in population-imbalanced case between a fermion with 
momentum k + q and spin f and a fermion with mo¬ 
mentum —k + q, and spin j. . As a result, the pair 
momentum is 2q. A finite pairing momentum implies a 
position-dependent phase of the order parameter, which 
in the Fulde-Ferrell (FF) case varies as a single plane 
wave A il)ia = Aqe l2q ’ ri i<5(r il -r i2 )M 

B. Mean field approximation 


The poles of the mean field single-particle Green’s 
function of the FF superfluidity in the present of SOC 
and the Zeeman field are defined by very long expres¬ 
sions. The numerical solution of the mean field set of 
two number equations, the gap equation and the equa¬ 
tion for the FF vector q is an ambitious task which will 
be left as a subject of our future research. 

In what follows, we consider only pairing between 
atoms with equal and opposite momenta, i.e. the BCS su- 
pcrfluid in the balanced system, and the Sarma^ super¬ 
fluid in the imbalanced case. In the imbalanced case, the 
Fourier transform of the zero-temperature single-particle 
Green’s function (1251) in the mean field approximation is 
given by2£ 


G m F (^j 


/ — £t(k) — h z — 2A (sinfca; + zsinfcj,) 0 

—2A (sin k x — i sin k y ) Ko m — £|(k) + h z — Ao 

0 — A 0 iu] m + ff(k) + h z 

Aq 0 — 2 A (sink x + isink y ) 


A 0 

0 

—2A (sinA; x — jsinfcy) 
JWm + £j.(k) - h z 


\ 


(26) 
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Here, the chemical potentials /r-pj. and the gap Ao are defined by the solutions of the mean-field number and gap 
equations: 


n = ft + fl = 1 - E 

E 


-im k)) 


x(k) 

fi(k) 


1 + 


5(k)+ry 2 (k) 


k L 


2 - /M k )) 


1 - 


V(5(k) + ^(k)) ( X 2 (k) + A 2 ) - 5(k)A 2 , 
S(k)+r] 2 {k) \ 


x(k) 

J w(k) r ^/(5(k) + rj 2 (k)) ( X 2 (k) + A 2 ) - 5(k)A 2 


( 27 ) 


nP = / t -A = -E 


1 


- /(«(k)) 


^(k) 

H(k) 


1 + 


X 2 (k) + Aq 


E 


k L 


2 - /M k )) 


1 - 


VWk) + V 2 (k)) (x 2 (k) + A 2 ) - 5(k)A 2 , 
X 2 (k) + A 2 \ 


^(k) 

J w (k) ^ V(5(k) + 7 ? 2 (k)) (y 2 (k) + A 2 ) - S(k)A 2 I ’ 


( 28 ) 


- = E 

u ^ 


k L 


2 - /(«(k)) 


1 


1 + 


7 ? 2 (k) 


+E 


k L 


2 - /M k )) 


J 2fi(k) ^ V(S(k) + r? 2 (k)) ( X 2 (k ) + A 2 ) - S(k)A 2 ^ 
1 /• 7 ? 2 (k) \ 


1 - 


2w(k) V ^(5(k) + ^(k)) ( x 2 (k) + Aq) - S(k)A§; ' 


( 29 ) 


Here, P = (/• j- — /j.)/(/f + /j,) is the polarization, /(w) is the Fermi-Dirac distribution function, and the following 
notations have been introduced: 

X(k) = \ (?t( k ) + £j.(k)), r?(k) = i (£ t (k) - £j.(k)) + h z , S( k) = 4A 2 (sin 2 k x + sin 2 fc y ) . 

The single-particle spectrum is 

H(k) =±\Js( k) + 1 (a ; 2 (k) + a ; 2 (k)) + 2^/(S(k) + ^(k)) ( X 2 (k) + A 2 ) - 5(k)A 2 , 


;(k) = ±\Js( k) + i (w2.(k) + w 2 (k)) -2^(S(k)+ rf (k)) ( X 2 (k) + Ag) - S (k) a|, 


where w±(k) = i/ X 2 (k) + Ag ± rj( k). 

In the mean-field approximation, the components G^ F 2 (ni,ri 2 = 1,2, 3,4) of the zero-temperature single-particle 
Green’s function GMF{k,id) are given by 


rtMF 

n-i n r> 


(k,w)= 


A ni n 2 (k) 


B n 


s (k) 


G„ in2 (k) 


+ 


D„ 


s (k) 


w — f 2 (k) + t0 + oj + fl(k) — *0+ id — w(k) T *0+ w+w(k) —*0+ 


( 30 ) 


where the corresponding expressions for A niV2 (k), B ni7l2 (k), C nin2 (k) and D niTl2 (k) can be obtained by inverting the 
matrix (l 26 l) . As an example, we have provided in the Appendix An(k),Bn(k),Gn(k) and Z?n(k). 


C. The Bethe-Salpeter equation for the collective excitations in the generalized random phase 

approximation 

The spectrum of the collective modes will be obtained by solving the BS equation in the GRPA. As we have 
already mentioned, the kernel of the BS equation is a sum of the direct Id = ST, F /6G and exchange I exc = 5T, H /SG 
interactions, written as derivatives of the Fock (l23l) and the Hartree (l24l) parts of the self-energy. Thus, in the GRPA 
the corresponding equation for the BS amplitude rli is given by 




rii n 3 
n 2 ri 4 


;(Q) 


Id 


U4 77-6 


^3 \ , T I n 3 n 5 

I A exc 


77-4 ^6 


n 6 ,n 5 5 


(31) 
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where I d ( Ul 713 = -ri 0 ) (ni,n 3 )A^ 0 irl 0 ) (n 4 ,n 2 ) and I exc ( ^ = iri 0 ) (?ii, n 2 )A^ 0 ]rl 0 ) (n 4 , n 3 ) are the di- 

\ 712 TI 4 J H K \ Tl 2 1 L 4 J H H 

rect and the exchange interactions, correspondingly. The two-particle propagator in the GRPA is defined as 
follows: 


A (0) 


n 1 n 3 
n 2 ri 4 


W (Q ) J — K nin3nin2 
Ani7i3 (k T Q) A n4n2 (k) 

w- (n(k + Q) + n(k)) 


Bn 


(k 


<in r d d k 

27r J (2ir) d 

Q) A, l4n2 (k) 


qMF 


(k + Q,fi + w(Q))G£* (k,fi) 


Cn\ri 3 (k + Q )D 

114712 (k) D nin 3 (k + Q)C'„ 4 n 2 (k) 


,(k + Q)£>n 4 n 2 (k) 


+ (S2(k + Q) - 
'n x n 3 (k + Q)A 


fi(k)) 

n 4 n 2 (k) 


oj — (w(k + Q) + ic(k)) oj T (u^(k T Q) -\- 

1)1] 7(3 (k T Q)A. 


Un x n 3 (k + Q)f3 n4Tl2 (k) 


w(k)) 

n 4 n 2 (k) 


(32) 


oj - (fl(k + Q) + w(k)) oj + (fi(k + Q) + w(k)) oj — (w(k + Q) + fl(k)) oj + (w(k + Q) + fl(k )) 1 

The BS equation ( 1311 ) can be written in the matrix form as ^I + UZ^j T = 0, where I is the unit matrix, the matrix 
Z is a 16 x 16 matrix, and the transposed matrix of T is given by: 


§ T = (T'p, 




Q 

1,2 




Q 

1,3 




Q 

1,4 


*£l 


'k 


Q 

2,2 


^ 2,3 




Q 

2,4 




Q 

3,1 


Vj/Q 

^ 3,2 




Q 

3,3 


*£4 




Q 

4,1 




Q 

4,2 




Q 

4,3 




The collective-mode dispersion, w(Q), follows from the condition that the secular determinant det\I + UZ\ is equal 
to zero. By standard algebraic manipulations, the 16 x 16 determinant, that follows from the BS equation, can be 
reduced to the following 12 x 12 secular determinant: 


^12x12 — 


1 + A /2 (A1221 - A1441) 

-AA1211 

AA1411 

-AA1121 

A /2 (Ann - A1331) 

AA1321 

A /2 (A2221 — A2441) 

1 - AA2211 

AA2411 

—AA'2121 

A /2 (A 2 iii — A2331) 

AA2321 

A /2 (A4221 - A4441) 

-AA4211 

1 + A A4411 

-AA4121 

A /2 (A4111 — A4331) 

AA4321 

A /2 (A1222 — A1442) 

— AA"i 212 

AA1412 

1 - AA1122 

A /2 (A m2 — A1332) 

AA1322 

A /2 (A1222 — A1442) 

— AA2212 

AA2412 

—AA2122 

1 + A /2 (A2112 — A2332) 

A A2322 

A /2 (A'3222 — A'3442) 

— AA3212 

A A'3412 

—AA3122 

A /2 (A3112 — A'3332) 

1 + A A3322 

A /2 (A'2223 — A'2443) 

— AA 2213 

A A'2413 

—AA2123 

A /2 (A'2113 — A'2333) 

A A2323 

A/2 (A'3223 — A3443) 

AA3213 

A A'3413 

—AA3123 

A /2 (A3113 — A3333) 

A A3323 

A /2 (A4223 — A4443) 

— AA4213 

AA4413 

—AA'4123 

A /2 (A4113 — A4333) 

A A4323 

A /2 (A1224 - A1444) 

-AAi 214 

AAi 414 

-AAu 24 

A /2 (A 1114 — A1334) 

AA1324 

A /2 (A3224 — A3444) 

-AA3214 

AA3414 

-AA3124 

A/2 (A 3 h 4 - A3334) 

A A3324 

A /2 (A4224 — A4444) 

— UK 4 2 i 4 

A A4414 

— AA'4124 

A /2 (A4114 — A4334) 

A A4324 


UK 1231 

UK 2 231 
AA 4231 
U K 12 32 
UK 2232 

UK 3232 

1 + U A 2233 

UK 3233 

U A'4233 
UK 1234 
UK 3234 
UK 4234 


U/2 (AA 441 - 
U /2 (A 24 41 - 
U /2 (A444I - 
u/2 (A u 42 - 
A /2 (A 24 42 - 
A /2 (A3442 - 
A /2 (A 24 43 - 
1 + U/2 (JT3443 
A /2 (A4443 - 
A /2 (Ad444 - 
A /2 (A3444 - 
A /2 (AT 4444 - 


AT1221) 
A2221) 
AI4221) 

A1222) 

A 2222) 

A 3222) 

A2223) 

— A3223) 
A4223) 
AI 1224 ) 
AI3224) 

A 4224) 


-AAT 143 i 
—AA'2431 
—A A/4431 

-aa 1432 

— AAT 2 432 
— AA 3432 
— AAT 2 433 
— AA 3433 
1 — A A 4433 
— AA 1434 

— AA'3434 

— UK 4434 


UKmi 
A Ad i 41 
UK4141 
A An 42 
A A 2142 
AA3142 
AA2143 

AA3143 

A A'4143 
1 + AATH44 
A A/3144 

AA4144 


1 - 


A AT1341 
AA 23 4 i 
A A4341 
A A, 342 
A A " 23 42 
A A3342 
A A2343 
A A3343 
A A4343 
AA1344 
AA3344 
A A'4344 


A /2 (ATi 33 i - 
A/2 (A'2331 — 
A/2 (A'4331 — 
A /2 (A/332 - 
A /2 (AT2332 — 
A /2 (AT3332 — 
A /2 (A2333 — 
A /2 (A3333 — 
A /2 (K 4333 - 
U /2 (A'i 334 — 
U /2 (A'3334 — 
1 + A /2 (A4334 


Ann) 

A2111) 

A 4111 ) 

A 1112 ) 

A2112) 

A 3112) 
A2113) 
A3113) 
A4113) 

A 1114 ) 

A3114) 

A 4114 ) 


(33) 


In the above secular determinant, there are 144 two-particle Green’s functions, but only 78 of them are different: 
Ann, A1144, A1122, A2211, A2222, A2233, A3322, A3333, A3344, A4433, A4411, A4444, and 66 different A niin2in3jra4 
{K nin2n3ni = K n 2 t n lt ni,n 3 )- It is worth mentioning that within the Gaussian approximatio n 30 : 31 the secular deter¬ 
minant includes only 8 of them, A'2233, A/3322, A1144, A4411, A'1234, A3421, A2413, and A1324: 


Z 2 x 2 — 


1 + -J (A2233 + A1144 — A1234 — A'2143) -J (A1414 + A'2323 — A2413 — A1324) 
y (A1414 + A 2 323 — A'4231 — A'3142) 1 + T (U -3322 + A'4411 — A ' 3 412 — A4321) 

I- 


( 34 ) 


In next Section, we shall see that the speeds of sound in the case of Sarma superfluid, calculated within the 
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Gaussian approximation and from the BS equation in the 
GRPA, are similar, but the dispersions around the roton 
minimum, provided by the two methods are remarkably 
different. 


III. NUMERICAL RESULTS 

A. Collective modes of an imbalanced 2D SOC 
Sarma superfluid in an optical lattice 

We use the secular determinants (l33ll and (l3dl) to cal¬ 
culate the collective-mode dispersion and the speed of 
sound of an imbalanced Fermi gas in a 2D square opti¬ 
cal lattice (lattice constant a) with a Rashba SOC. We 
consider a weak coupling limit U = 2.64J, where the BS 
equation in the GRPA provides a good approximation for 
the collective-mode energies. The strength of the SOC is 
A = 0.1 J, and the system parameters are chosen so that 
the density of the majority and minority components are 
/j- = 0.275 and /j, = 0.225, respectively (the correspond¬ 
ing polarization is P = 0.1). The chemical potentials 
/r-|- = 2.857J, = 2.186J and the gap A = 0.266J are 
obtained by numerically solving the number and the gap 
equations (Ha-®. 

In Fig. 0 we plot the collective-mode excitation en¬ 
ergy ui(Q x ) as a function of the wave vector Q = ( Q x , 0) 
along the x-axis, calculated by the Gaussian approxima¬ 
tion (triangles) and by the BS approach in the GRPA 
(circles). The speed of sound, provided by the Gaus¬ 
sian approximation, is u = 1.66 Ja/h , while the speed 
of sound, calculated within the BS approach is u = 1.35 
Ja/H. Thus, the Gaussian approximation overestimates 
the speed of sound by about 23%. More interestingly, the 
dispersion curve calculated from the BS equation clearly 
shows the existence of a roton-like minimum, while there 
is no such a minimum within the Gaussian approxima¬ 
tion. The corresponding roton gap is A r = 0.2025J and 
the critical flow velocity, obtained around the roton min¬ 
imum from the BS equation, is v c = 0.51 Ja/h. As in 
2D case, shown in Fig. [2] the two dispersion curves are 
remarkably different around the roton minimum; instead 
of the expected roton-like structure, the dispersion curve 
provided by the Gaussian approximation monotonically 
increases in this region. At higher wave vectors, Q x > 
0.16 7r/a, the two dispersion curves have essentially the 
same behavior. 

Our results suggest that the Gaussian approximation 
overestimates the speed of sound of the Goldstone mode, 
and fails to reproduced the roton-like structure of the 
collective-mode dispersion which appears after the linear 
part of the dispersion. The question naturally arises here, 
whether the Gaussian approximation still can be used to 
estimate the speed of sound if one takes into account not 
only the SOC, but Zeeman fields as well. This question 
will be answered in the next subsection, where we inves¬ 
tigate the speed of sound of a balanced superfluid Fermi 
gas in a two-dimensional square optical lattice with a 



FIG. 3: Collective modes dispersion oj(Q x ) of an imbalanced 
SOC Fermi gas in a 2D square optical lattice along the (Q x , 0) 
direction, obtained by the Gaussian approximation (triangles) 
and the Bethe-Salpeter equation in the GRPA (circles). 



h z - VA 2 +|i 2 (units of J) 


FIG. 4: The chemical potential (diamonds) and the gap 
(squares) of a balanced superfluid Fermi gas in a 2D square 
optical lattice as a function of the out-of-plane Zeeman field. 
The system parameters are: filling factor / = 0.5, the attrac¬ 
tive interaction U = 5.2 J, and the strength of the Rashba 
spin-orbit-coupling A = 0.1J. 


Rashba spin-orbit coupling (in the xy plane) and an out- 
of-plane Zeeman field. 
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h z - VA 2 +|T (units of/) 


FIG. 5: The speed of sound along the x direction as a function 
of the Zeeman field, calculated within the Gaussian approx¬ 
imation (circles) and by the BS approach (triangles). The 
values of the chemical potential and gap are shown in Fig. [3] 


B. Speed of sound near the transition from the 

gapped superfiuid phase to the topological phase 

We consider a system with a filling factor / = 0.5, 
and an attractive interaction U = 5.2 J. The strength 
of the Rashba spin-orbit-coupling is A = 0.1 J, and there 
exists an out-of-plane Zeeman field h z . In such a system 
a phase transition can be accessed by varying the Zeeman 
field. The transition from the gapped superfluid phase to 
the topological phase is characterized by the quasiparticle 
excitation gap that closes at h c = i//i 2 + A 2 and reopens 
with increasing h z > h c . 

In Fig. [4] we present the chemical potential (dia¬ 
monds) and the gap (squares) for different Zeeman fields, 
obtained by solving numerically the number and the gap 
equations. The critical field h c marks the phase transi¬ 
tion between the topologically nontrivial (negative side) 
and the topologically trivial (positive side) superfluid 
phases. As can be seen, during this transition, the gap A 
is still finite even though the quasiparticle excitation gap 
is closed. This suggests that there is a quantum phase 
transition separating the parameter regimes h z < h c and 
h z > h c , even though the system in both regimes is an 
s-wave superfluid. 

In Fig. 0 we have plotted the speed of sound along x 
direction (Q = (Q x , 0)) as a function of the Zeeman field, 


calculated within the Gaussian approximation and from 
the BS equation in the GRPA. As can be seen, close to the 
phase transition boundary the speed of sound calculated 
within the two approaches is essentially the same. The 
inset in the figure shows that the minimum of the speed 
of sound is located at the phase transition boundary h c . 
The same behavior was previously found by applying the 
Gaussian approximation in the case of a 2D superfluid 
atomic Fermi gas with Rashba-type spin-orbit coupling 
and an out-of-plane Zeeman field) 28 : 30 and in the case 
of a 3D FF type of superfluid Fermi gas with Rashba 
spin-orbit coupling (in the xy plane) and two Zeeman 
fields [in-plane (hx) and out-of-plane (hz)]«21 Thus, our 
calculations are in agreement with the suggestion made 
in Ref. [31 | , that by measuring the minimum of the speed 
of sound one can unambiguously detect the topological 
phase transition boundary. 


IV. SUMMARY 

In summary, we have derived the BS equation in the 
GRPA for the collective excitation energy of a Fermi gas 
in a 2D square lattice with an attractive contact inter¬ 
action, assuming the existence of a nonvanishing Rashba 
SOC and an out-of-plane Zeeman field. We have calcu¬ 
lated the collective-mode dispersion within the Gaussian 
approximation, and from the BS equation assuming the 
existence of a Sarma superfiuid state. It is found that the 
Gaussian approximation: (i) overestimates the speed of 
sound of the Goldstone mode, and (ii) fails to reproduce 
the roton-like structure of the collective-mode dispersion 
which appears in the BS approach. 

We also have investigated the speed of sound of a 
balanced spin-orbit-coupled atomic Fermi gas near the 
boundary of the topological phase transition driven by 
an out-of-plane Zeeman field. It is shown that the mini¬ 
mum of the speed of sound is located at the topological 
phase transition boundary, and this fact can be used to 
confirm the existence of a topological phase transition. 
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APPENDIX 

The symbols I a ,b and J Qi t> at nonzero temperatures are 
defined as: 
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1 * - 

Ia ' b = 2 ak Q ° k - 


Q 


^ a,b ~ 2N ak -Q^ k -Q 


1 - f (a>_(k)) - / (a> + (k + Q)) 
w + f2(k, Q) - e(k, Q)] 

■l-/(w_(k))-/(w+(k + Q)) 


■n(k,Q)-e(k,Q)] 


l-/(a; + (k))-/(u;_(k+Q)) ' 
ui + fi(k, Q) + e(k, Q)] 

1 ~/(^+( k )) ~/(^-(k + Q)) ~ 
w + H(k, Q) + e(k, Q)] 


Here, e(k, Q) = £7(k + Q) + E{ k), H(k, Q) = r?(k) - ^(k + Q), E( k) = Vx 2 ( k ) + A 2 , w±(k) = £(k) ± r?(k), and a 
and 6 are one of the following form factors: 

7k,Q = Wk«k+Q + Ck'Ck+Q, ^k,Q = Wk'«k+Q — t’kt’k+Q, 7k,Q = Mk^k+Q — Wk+Cp’k, Wk,Q = Uk'^k+Q + Wk+Ql’k, 



The functions A nin2 (k), B nin2 (k), C ni n 2 (k) and D nira 2 (k) are easily obtained by inverting the matrix (l26l) . Here, 
we provide the expressions for An(k), Bu(k), Cn(k) and Z?n(k) (the k-dependence of ^^(k), 5 (k), H(k), and w(k) 
is not explicitly shown): 


An(k) = 


[~h z + a(k)] - [Ao + S( k) + [h z - &(k)] [h z + £ t (k)]] - [A 2 + 5(k) + [h x - &(k)]] + [h z + £ t (k)] H 2 (k) + H 3 (k) 

2 f 2 (k) [ft(k) - w(k)] [fi(k) + w(k)] 


B n(k) = 


l-h z + a(k)] - [Ao + S( k) + [h z - &(k)] [h z + £ t (k)]] - [A 2 + S(k) + [h z - &(k)]] Q — [h z + £ t (k)] H 2 (k) + H 3 (k) 

2 H(k) [H(k) - w(k)] [H(k) + w(k)] 


C'n(k) = 


[~h z + £+(k)] - [A 0 + S(k) + [h z - ^(k)] [h z + £ t ( k )]] + [A 2 + S( k) + [h z - ^(k)]] u - [h z + g t (k)] u; 2 (k) - u; 3 (k) 

2w(k) [fi(k) — w(k)] [H(k) + w(k)] 


D n(k) = 


[~h z + gj,(k)] - [A 0 + S(k) + [h z - ^(k)] [h z + g t (k)]] + [A 2 + S( k) + [h~ - ^(k)]] u + [h~ + g t (k)] u; 2 (k) - u; 3 (k) 

2w(k) [fi(k) — w(k)] [H(k) + w(k)] 
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